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ABSTRACT 



We perform global unstratified 3D magnetohydrodynamic simulations of an astro- 
physical boundary layer (BL) - an interface region between an accretion disk and a 
weakly magnetized accreting object such as a white dwarf - with the goal of under- 
standing the effects of magnetic field on the BL. We use cylindrical coordinates with 
an isothermal equation of state and investigate a number of initial field geometries in- 
cluding toroidal, vertical, and vertical with zero net flux. Our initial setup consists of 
a Keplerian disk attached to a non-rotating star. In a previous work, we found that 
in hydrodynamical simulations, sound waves excited by shear in the BL were able to 
efficiently transport angular momentum and drive mass accretion onto the star. Here 
we confirm that in MHD simulations, waves serve as an efficient means of angular mo- 
mentum transport in the vicinity of the BL, despite the magnetorotational instability 
(MRI) operating in the disk. In particular, the angular momentum current due to waves 
is at times larger than the angular momentum current due to MRI. Our results suggest 
that angular momentum transport in the BL and its vicinity is a global phenomenon 
occurring through dissipation of waves and shocks. This point of view is quite different 
from the standard picture of transport by a local anomalous turbulent viscosity. In 
addition to angular momentum transport, we also study magnetic field amplification 
within the BL. We find that the field is indeed amplified in the BL, but only by a factor 
of a few and remains subthermal. 

Subject headings: accretion, accretion disks - magnetohydrodynamics — waves - insta- 
bilities 



1. Introduction 



The mechanism of angular momentum transport in the boundary layers (BLs) of accreting 
compact objects and stars is an open problem that is both interesting from a physical point of view 



and important for modeling BL spectra from these sources (Narayan 1992; Popham &; Narayan 



1992, 1995; Inogamov fc Sunyaev||1999 ; Piro &; Bildsten 2004). In this work, we shall focus on the 
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case of an accretion disk that does not undergo disruption by the magnetic field of an accreting 



opbject (for the case of magnetically disrupted disk see Ghosh &; Lamb (1978), Koldoba et al. 
( |2002| )) and extends all the way to the surface of the compact object or star. Such a situation 
is expected to be naturally realized in weakly magnetized neutron stars producing Type-I X-ray 
bursts as well as during FU Orioni outbursts in pre-main sequence stars and dwarf nova outbursts 
in cataclysmic variables when the mass accretion rate is high. 

Accretion of matter through the BL must be accompanied by the angular momentum transport 
in the layer, but the nature of this transport has remained elusive. The BL is linearly stable to 



the magnetorotational instability (MRI; Velikhov (1959); Chandrasekhar (1960); Balbus &; Hawley 
( |1991D ) that is thought to give rise to anomalous viscosity in the bulk of well-ionized accretion disks, 
because the angular velocity, fi, necessarily increases with radius in this region. Of course linear 
stability does not imply a nonturbulent flow. Turbulence, once generated by some nonlinear mech- 
anism, can persist without decaying in certain classes of rotating, linearly stable hydrodynamical 
flows ( Lesur fc Longaretti|2005| ) . However transport mechanisms, both magnetic and hydrodynamic 
involving an anomalous viscosity have never been definitively demonstrated to operate inside the 
BL. 



Recently Belyaev &; Rafikov (2012) and Belyaev et al. ( 2012a|b ) showed both theoretically and 
computationally that non-axisymmetric acoustic modes are excited in the BL when it is thin in 
comparison with the stellar radius. These acoustic modes manifest themselves as sound waves that 
propagate away from the boundary layer, efficiently transporting angular momentum into both the 



star and the disk and driving accretion. However, the studies of Belyaev & Rafikov (2012), Belyaev 
eTaTI ( [2012a|b| ) were purely hydrodynamical, and it is the purpose of the present study to check 
their results on angular momentum transport by waves generated in the BL in the presence of a 
magnetic field. 

There are several reasons why the MHD case could be fundamentally different from the hy- 
drodynamical case. First, the presence of a seed magnetic field in the disk leads to turbulence 
generated by the MRI. Potentially, the turbulence could interact with the acoustic modes in the 
disk, modifying their properties or washing them out completely. Second, magnetic field advected 



into the BL by accreting gas will be amplified by the shear present there (Pringle 1989; Armitage 



2002). A strongly magnetized BL could have different properties than an unmagnetized one, again 



influencing the properties of acoustic modes. Third, magnetic field in the BL can give rise to 
other MHD effects such as the Tayler-Spruit dynamo ( |Tayler|[T973| |Spruit|[l999| |2002[ ), which may 



produce additional momentum transport in the BL (Piro Sz Bildsten 2007). 



One of our aims in this work is to determine the extent to which a magnetic field affects the 
properties of and angular momentum transport by waves generated in the BL. We find that in the 
presence of a magnetic field, the acoustic modes studied by Belyaev et al. ( 2012a|b ) become magne- 
tosonic modes. However, when the ratio of gas pressure to magnetic pressure f3 ^> 1, the magnetic 
field has only a weak, O effect on the dynamics of acoustic modes. In our simulations, (3 is 
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indeed large, and we find that magnet osonic modes play an important role in angular momentum 
transport in the BL and its vicinity. 



Previously, Armitage (2002) and Steinacker & Papaloizou (2002) each performed MHD simula- 



tions of the BL. However, these studies did not reveal the importance of the transport by waves, and 



we discuss our work in relation to theirs in §8l Also, Heinemann & Papaloizou (2009a), Heinemann 



& Papaloizou (2009b), Heinemann 8z Papaloizou (2012) studied angular momentum transport by 



spiral density waves in MRI turbulent disks, which were sourced by vorticity generated in non- 
linear MRI turbulence. However, the mechanism for driving these waves is entirely different from 
the mechanism for generating the magnetosonic modes discussed here. The latter are sourced by 
shear in the BL and launched away from it into both the star and the disk. 

There have also been a number of other analytical and numerical studies of disk wave phe- 
nomena in the vicinity of the innermost stable circular orbit of a black hole and near the BL (Li 
et al.||2003[[Li fc Narayan||2004[ [LaTfc Tsang||2009{ |Tsang feT^|2009a|b|c| |Fu fc Lai||2012| ). How- 
ever, these studies were primarily addressing the possible connection between the non-axisymmetric 
modes and quasi-periodic oscillations observed in accreting objects, whereas our main focus is on 
the fundamental nature of the angular momentum transport within the BL. Moreover, with the ex- 



ception of Tsang &; Lai (2009b) and Li & Narayan (2004), these studies did not include a magnetic 
field. 

The paper is organized as follows. In ^2j we discuss the technical aspects of our numerical 
model, and in ^3] we introduce notation and define quantities used throughout the paper when 
discussing angular momentum transport. Then, in ^4] we talk about angular momentum transport 
by MRI in the disk and show that our results agree with previous studies of the MRI. £[6] discusses 
how acoustic modes studied by Belyaev et al.| ( |2012a|b ) are modified by the presence of a magnetic 
field, becoming magnetosonic modes. Finally, §[7] discusses angular momentum transport by these 
magnetosonic modes in the star, BL, and inner disk. 



Numerical Model 



We begin by summarizing the equations and the numerical model used in our simulations. In 



this paper, we use a numerical model that is identical to the one presented in |Belyaev et al. (2012b), 
except for the addition of a magnetic field. Thus, we will only briefly summarize the aspects of 



the numerical model that do not relate to the magnetic field, referring the reader to Belyaev et al 



(2012b) for rest of the details. 
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2.1. Governing Equations and Nondimensional Quantities 



For our simulations, we use Athena (Stone et al 2008) to solve the MHD equations with an 



isothermal equation of state in cylindrical geometry. These equations are 

^ + V-(pv) = (1) 
9W +V^(H = -V (p+^) + -(B- V)B-pV$ (2) 



dt \ 2[i J [i 

— = Vx(vxB) (3) 
P = ps 2 . (4) 

Here, v is the velocity, p is the density, P is the pressure, $ is the gravitational potential, s is 
the sound speed, and B is the magnetic field. We also use w throughout this work to denote the 
cylindrical radius. 

Adoption of an isothermal equation of state is equivalent to assuming fast cooling to a set 
temperature. Although this is not a realistic assumption, it allows us to simply attain a quasi- 
steady state and focus on exploring the effect of a magnetic field on the dynamics of the BL. 

We nondimensionalize quantities by setting the radius of the star to vj* — 1, the Keplerian 
angular velocity in the disk just above the surface of the star to Q^^m*) = 1, the density in the 
disk, initially a constant, to p = 1, and the magnetic permeability to ji = 1. Finally, we define a 
characteristic Mach number for our simulations as M = 1/s. This corresponds to the true Mach 
number of fluid rotating at the Keplerian velocity at w — w*. 



2.2. Summary of Aspects of the Numerical Model not Relating to the Magnetic 

Field 

The simulation setup consists of a radially stratified, unrotating star attached to a Keplerian 
disk via a thin interfacial region. The simulations are performed in cylindrical geometry with 
an isothermal equation of state and are not vertically stratified. The interfacial region smoothly 
connects the rotation profile of the star to that of the disk and is chosen to be as thin as possible, 
while still being resolved with ~ 10 cells. The highest density in the star is approximately 4 x 10 6 
the density in the disk. Thus, the inner regions of the star have a high enough inertia to avoid 
being spun up over the course of a simulation. 

We do not incorporate self-gravity into our runs, but instead use a fixed gravitational potential 
V — —1/w. We initialize each of our simulations to be in hydrostatic equilibrium and add random 
perturbations to the radial velocity from which instabilities develop. For additional information on 



aspects of our numerical model that do not pertain to the magnetic field see Belyaev et al. (2012b). 
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2.3. Initial Magnetic Field Configuration 

We now discuss the initial magnetic field strengths and geometries used in our simulations. 
For all our simulations, we set B — for w < (an unmagnetized star), and initialize one of 
three possible magnetic field geometries for w > w*\ 

%Lz, NVF 
B(m>w in z) = { NAF . ( 5 ) 



^sin 



2tt 



{w — w+ 



ZNF 



NVF simulation have a constant vertical magnetic flux through any radius, NAF simulations have 
a constant azimuthal magnetic flux, and ZNF simulations have on average zero net flux. For the 
ZNF simulation, we set \ w — Az/2, where Az is the extent of the simulation domain in the 
vertical direction. Thus, we have a rapid variation of the field in the radial direction, which allows 
comparison with shearing box ZNF simulations. 

In all of our runs, we use an initial value of Bq = .002 in units for which the magnetic 
permeability is (i = 1. Defining the plasma parameter /3 as 

the initial value of (3 using Bq for a baseline is 0q ~ 6200. Note that in the NVF and ZNF 
simulations, (3 varies as a function of position and fio specifies the minimum initial value of f3 in 
the disk. 

The particular value of Bq = .002 is motivated by two considerations. The first is the need to 
resolve the wavelength of the fastest growing mode of the MRI in the case of a pure vertical field 
(NVF simulations), and the second is the need to fit several wavelengths of the fastest growing 
mode within the simulation domain. 

For a global disk simulation in which the domain has constant width and resolution in the 
z direction, these two criteria are generally at odds with one another. This is because the z- 
wavelength of the fastest growing axisymmetric mode for constant vertical field in a Keplerian disk 



is given by Balbus & Hawley (1991) 



Wi-Vl^- (7) 

For a Keplerian disk with constant B z field and constant p, A^mri uj 3 / 2 . For our NVF field 
geometry, the magnetic field strength falls off as 1 jvo so A^mri w 1 / 2 . This is a weaker dependence 
than for a constant field and is the main reason we use a constant flux criterion in the NVF case, 
rather than simply using a constant field. 

For the ZNF case, we no longer have to worry about fitting the vertical wavelength of the 
fastest growing mode within the simulation domain, since the field strength varies as a function of 
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radius and the MRI grows on all scales. The same goes for NAF simulations, because instability 
manifests itself differently if the initial field is azimuthal rather than vertical. In the former case, 
the most unstable modes are non-axisymmetric ones that undergo transient exponential growth in 
the shearing sheet approximation ( Balbus fc Hawley||1992 ; Tagger et al.| [T992) or pure exponential 
growth in a global disk analysis ( Matsumoto fc Taj ima| 1995 ; Coleman et al.|1995 ; Ogilvie &; Pringle 



1996; Terquem & Papaloizou 1996). In a shearing sheet analysis, the modes which undergo the 



most vigorous growth are the ones that have the largest values of k z for a given value of k x and k y 
(Balbus &; Hawley 1992; Hawley et al. 1995). Thus, resolution of these modes in a simulation is 
limited by the grid resolution, but their vertical wavelength is always able to fit within the vertical 
extent of the simulation domain. 



2.4. Simulation Parameters 



The only physical parameter varied across simulations in this study is the initial field geome- 
try. In principle, other parameters such as the Mach number, the vertical extent of the simulation 
domain, the presence of vertical stratification, etc., would also have an important effect on the 
dynamics of the BL. However, it is computationally expensive to undertake an exhaustive parame- 
ter study, so we limit ourselves to varying only the initial field geometry in this exploratory work. 
This allows us to explore the major effects a magnetic field has on BL dynamics with no added 
frills. We mention that Belyaev et al.| ( |2012a|b| ) have already explored how varying the dimen- 
sionality (2D/3D), stratification (stratifled/unstratifled), and Mach number affects the outcome of 
hydrodynamic BL simulations. 

Some important parameters that are fixed across different runs are the Mach number, set to 
M — 9; the vertical extent of the simulation domain, set to Az = s/Qk(j&*) = 1/9; the (j) extent 
of the simulation domain, set to Acp = 2tt/7; and the (p resolution set to = 384 cells. Table [l] 
lists important numerical parameters that vary from simulation to simulation. 

The specific choice of = 27r/7, represents a compromise between computational resources 
and simulating a wedge that spans a wide enough azimuthal angle. Also, setting = 27r/7, 
we resolve all modes having m a multiple of seven. The m = 14 mode happens to be close to a 

9, which means its 
14 mode 



geometric resonance (equations (44)-(45) of Belyaev et al. (2012b)) for M 
amplitude is enhanced relative to modes with a similar m number. Therefore, the m 
typically peaks in our simulations, simplifying the analysis, since we can study a single dominant 
mode rather than a superposition of modes having similar amplitudes. Note that the finite extent of 
our simulations in the azimuthal direction precludes us from exploring the Tayler-Spruit dynamo 



(Tayler 1973; Spruit 1999, 2002) in the BL, which would require long azimuthal wavelength to 



be resolved. However, the very existence in the BL of this dynamo previously shown to operate 
only in the incompressible limit should be seriously questioned in a highly dynamic, compressible 
environment such as the BL. 
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The boundary conditions we use are periodic in the z and <j) directions, and do-nothing at the 
inner radial boundary. A do-nothing boundary condition means that the fluid quantities retain 
their initial values on the boundary for the duration of the simulation. This is preferred over an 
open (outflow) boundary condition for the inner boundary, since it ensures that the star does not 
slowly "drift" out of the simulation domain. It is also preferred to a reflecting boundary condition, 
since it is (partially) transparent to outgoing waves ( |Belyaev et al.|[2012b ). 

At the outer radial boundary, we implement two different types of boundary conditions. The 
first type, which we call OBCO, is simply an open (outflow) boundary condition. The second type, 
OBC1, is a do-nothing boundary condition, just like at the inner radial boundary. In simulations 
implementing an outer do-nothing boundary condition, we initialize a nonmagnetic "buffer zone" 
between a fiducial radius in the disk n7B, m ax and the outer boundary. The magnetic field is initialized 
to be zero in the region zu^q m ax 5^ ^ 5^ ^outj where tu ou t denotes the radius of the outer boundary. 
In all of our simulations implementing OBC1, n7B,max ~ ^. 



1 
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Fig. 1. — Panels a,b show the magnitude of the magnetic field in simulation M9d at times t = 
and t = 500, respectively. 

To demonstrate the buffer zone more clearly, we plot the magnetic field magnitude from simula- 
tion M9d at t = and t = 500 in panels a and b of Fig. [I] In panel a, the field decays oal/w between 
the star at w — 1 and the buffer zone at vo = ti7B,max = 4. In the region zz7B, m ax < w < ^out, the 
field is initially set to zero. In panel b, magnetic field has penetrated into the buffer zone, as a result 
of the natural evolution of the MHD equations, but has not yet reached the outer boundary. In 
general, simulations implementing OBC1 are stopped before the magnetic field can diffuse through 
the buffer zone and reach the outer boundary. This allows us to run simulations for a long time, 
without having to worry about the effect of magnetic field at the outer boundary on the rest of the 
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simulation domain. 



label 


B type 


outer BC 


ro-range 


N w 


N z 


^run 


M9a 


NVF 


OBCO 


(.85, 6.1) 


4096 


128 


450 


M9b 


NAF 


OBCO 


(.85, 7.85) 


6144 


64 


450 


M9c 


ZNF 


OBCO 


(.85, 6.1) 


4096 


128 


450 


M9d 


NVF 


OBC1 


(.85, 6.1) 


4096 


64 


800 


M9e 


NAF 


OBC1 


(.85, 6.1) 


4096 


48 


1000 


M9f 


ZNF 


OBC1 


(.85, 6.1) 


4096 


64 


1800 



Table 1: Summary of simulation parameters for MHD simulations. The columns from left to right 
are: simulation label, initial field geometry (net vertical flux, net azimuthal flux, or zero net flux), 
outer boundary condition, radial extent of the simulation domain, number of grid points in the 
radial direction, number of grid points in the vertical direction, and the total time for which the 
simulation was run. 



3. Angular Momentum Transport: Equations 

In this section, we introduce the notation used throughout this paper when discussing angular 
momentum transport. When analyzing simulation results, we perform both density weighted and 
volume weighted averages of quantities over the z and (f) directions. A volume weighted average is 
denoted by an overline: 

where Az and Acf) are the extents of the simulation domain in the z and <j) dimensions respectively. 
A density weighted average is denoted by angle brackets: 

(/) = j— \ dzfypf, (9) 



where = p is the volume weighted average density. 

For ideal MHD, the one dimensional equation of angular momentum transport can be formu- 
lated in terms of density weighted and volume weighted averages as ( |Balbus fc Papaloizou|[l999[ ) 



-(w 2 j; n) + -—(w 3 nj; {v w ) + w 2 T^) = o, (10) 



where £1 = (v^/w is the density -weighted angular velocity, and T^ds is the nj, (j) component of the 
stress tensor 

T m( t> = PVzu - ^) - BvjB() . (11) 

p 
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The first term in the space derivative in equation (10) equals the advected angular momentum 



current, C% divided by 27r, and the second term equals the stress angular momentum current, C5, 
divided by 2tv. For a steady state disk, Ca + Cs — 0, although the disks in our simulations are not 
in steady state, and undergo temporal evolution. 

The stress current, C#, can be further decomposed into hydrodynamical and magnetic terms 

Cs^fat) = 27rw 2 ^ (v m (v^ - wQ)) (12) 



-2ttw z 



(13) 



C s = Cs,h + Cs,b. (14) 
Note that Cs includes contributions to the stress from both waves and the action of turbulence. 

Finally, we define an a parameter to parametrize the angular momentum transport rate by 
turbulent viscosity. There are many possible definitions for a, especially in the context of the BL 
( Narayan||1992 ; Popham fc Narayan||1995| ), but the ones used in this work are 



a(zu, t) = 



a v {w,t) 



(15) 



(16) 



TiQsHdQ/dlnvci 

The first form of a is just the ratio of the stress to the pressure. The second is the a for a 



turbulent anomalous viscosity (Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974), where H 



is the vertical scale of the problem. In an astrophysical disk, H = s/ft and increases radially for 
constant s. Our simulations, however, are in cylindrical geometry and H = Az, where Az, the 
extent of the simulation domain in the vertical direction, is a constant. 



Analogous to the splitting of Cs into hydrodynamical and magnetic terms (equation [14]), it 
is useful to split a into hydrodynamical and magnetic components: 

(Vm (V<l> ~ G7f2)) 



a H (m,t) = 



a — an -\- 



(17) 

(18) 
(19) 



Although a is a useful parametrization of angular momentum transport by MRI turbulence 
in the disk, it is not applicable to transport by waves. In the latter case, the relevant quantity is 



Cs (Belyaev et al. 2012b), and since angular momentum is transported by both waves and MRI 



turbulence in our simulations, we need to consider both Cs and a when analyzing results. 

Another useful parameter we refer to often when analyzing simulations is /3 (equation [6]), and 
from now on we take 



/3(07,t) = =^ 



En* 2 



B 2 /2n 



(20) 
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4. Angular Momentum Transport by MRI in the Disk 

In this section, we focus on angular momentum transport in the disk and show that we obtain 
resolved MRI turbulence in our simulations. Thus, this section serves as a check of our simula- 
tions against previous work, and the reader who is primarily interested in new results about wave 
transport of angular momentum may choose to skip ahead to the next section. 



4.1. Magnetic Pitch Angle 



The primary diagnostic we use to check that the disk reaches a resolved quasi-steady state of 
MRI turbulence is the magnetic pitch angle 



e B = 



sin 



- 1 (olbP) 



(21) 



Equation (21) gives an approximation for the angle of the magnetic field with respect to the (j) 
direction if B& > B w and B w > B z (IGuan et al.||2009|). ISorathia et all (120121) found that 9 B 



was a particularly good indicator of convergence, since it varies monotonically with resolution and 
converges to a value of 9 b ~ 13° in resolved simulations. A similar indicator of convergence is 
the relation a/3 « 1/2 (Guan et al. 2009; Hawley et al. 2011; Sorathia et al. 2012). However, a 



includes both hydrodynamical and magnetic stresses, whereas olb is calculated using only magnetic 
stresses. We find in our simulations that the hydrodynamical stress due to waves is non-negligible 
in the disk, but waves do not contribute to the magnetic stress. Thus, a diagnostic based on ^ 
is strongly preferred to one based on a, when assessing the convergence of MRI turbulence in a 
simulation that includes a BL. 

Fig. [2] shows radius time plots of 9b for simulations M9d,e,f in the annulus 1 < vj < 3. 
Ignoring the region near w — — 1 where the rotation profile is strongly non-Keplerian, we see 
that 9b ~ 12 — 14° in the disk proper for the NAF and ZNF simulations (panels b and c) once 
MRI turbulence is fully developed. In the NVF simulations, 9b ~ 12 — 14° in the inner disk, but 
rises with radius. Nevertheless, the values of 9b observed in our simulations are in good agreement 



with Sorathia et al. (2012). This gives us confidence that we have properly resolved the MRI in 
the disk. 



4.2. Convergence Across Different Resolutions 

The second diagnostic of convergence we perform is a measurement of the r.m.s. volume 
weighted average of B w ^ £>^, and B z for simulations with the same initial magnetic field geometry, 
but having a different resolution. Panels a,b,c of Fig. [3] show r.m.s. averages of B m , B^, and B z 
at t = 400 for NVF simulations (M9a,d), NAF simulations (M9b,e), and ZNF simulations (M9c,f), 
respectively. We see that the NVF and NAF simulations are converged, although the magnetic 



- 11 - 




Fig. 2.— Panels a,b,c show 6 B for simulations M9d (NVF), M9e (NAF), and M9f (ZNF). The band 
on the left hand side of panel c is due to a transient burst of activity that occurs when the MRI 
first goes nonlinear. Note the different scale on the time axis in the figures. 

field in the outer disk of the NAF simulation is still approaching a steady state. The reason why 
we did not perform the convergence study at a later time is that it was too expensive to run the 
high resolution simulations beyond a time of t — 450. 

In the ZNF case, magnetic field components for the low resolution simulation have a higher 
amplitude than for the high resolution simulation, and we have checked that this is true at all 
times, not just at t — 400. However, in (j |4.1| ) we showed that MRI turbulence is resolved in ZNF 
simulations based on the pitch angle indicator. 

Taken together, these two pieces of information suggest that MRI turbulence in our ZNF 
simulations is resolved, but that for a given simulation, the r.m.s field strength in the saturated 
state of MRI turbulence is influenced by the effective numerical viscosity and resistivity, which 
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Fig. 3. — The volume weighted r.m.s magnetic field components at t = 400 for NVF simulations 
(panel a), NAF simulations (panel b), and ZNF simulations (panel c). In each case, the solid line 
is the simulation with the higher z-resolution. 



depend on the grid resolution ( Fromang fc Papaloizou|[2007 ) . In order to make the results of ZNF 
simulations independent of resolution, we could add explicit viscosity and resistivity (Fromang et 



al. 2007). 



However, given the exploratory nature of this paper, the exact value of the viscous and resistive 
coefficients is of secondary importance. What is more important is that for a given value of these 
coefficients, MRI turbulence is properly resolved, which is indeed the case according to the pitch 
angle indicator of §4.1| 



4.3. The Value of a due to MRI Turbulence in the Disk 



The most important parameter that characterizes transport by MRI turbulence in the disk is 
a. However, a in our simulations is not a constant and is a function of both time and radius. One 
reason for this is that the timescale for instability is longer in the outer parts of the disk than in the 
inner parts, so the inner parts reach a quasi-steady state faster. However, even in the quasi-steady 
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state, we find a oc cu -3 / 2 , which is not a surprising scaling, and implies a constant value of a u 
(equation [l6]). The reason a and ol v have a different radial dependence in simulations is that the 
height of the simulation domain, Az is fixed, whereas the scale height in a Keplerian disk goes as 
H = s/Q oc iu 3 / 2 for constant s. Thus, we define 

a eff = aw 3/2 . (22) 

Since, s/Q = Az at w — — 1, a e R can be thought of as the "effective" value of a for a Keplerian 
disk, based on our simulation results. The assumption here is that the angular momentum transport 
rate is proportional to the vertical scale of the problem, which is fixed for our simulations but goes 
as H oc iu 3 / 2 for a Keplerian disk with a constant value of s. 

Another complication to measuring the a due to MRI turbulence is the presence of waves 
emitted from the BL, which contribute to the value of a in the disk. Therefore, rather than 
showing a e ff directly, we show o^effj which is simply the component of a e & that is due to magnetic 



stresses (i.e. a^eff — c^b^ 3 / 2 ). As mentioned in £4.1, «b is insensitive to waves, which contribute 



primarily to the hydrodynamical stress. Studies of the MRI typically find that a w 5/4a^ (Sorathia 



et al. 2012| ), which also roughly holds for our simulations. 



Fig. [4] shows radius time plots of c*B,eff f° r simulations M9d,e,f. Once a quasi-steady state has 
been reached, the value of «B,eff f° r the NVF simulation (M9d) is c*B,eff ~ -02 — .03, for the NAF 
simulation (M9e) is a B , e ff ~ -0075 - .0125, and for the ZNF simulation (M9f) is a B ,eff ~ -004 - .01. 



These values of a are consistent with those obtained by |Steinacker fc Papaloizou| ([2002]) in their 



BL simulations. However, the value of a in ZNF simulations does depend on the resolution, for the 
same reasons as the r.m.s. magnetic field value depends on the resolution (j|4.2|). 



5. Magnetic Field Amplification in the BL 

It has been speculated that the strong shear present in the BL can amplify the magnetic field 
to the point that f3 < 1 ( Pringle||l989 ). Armitage (2002) observed amplification of the field in the 



BL in MHD simulations with an isothermal equation of state, but found < .1 — .2, everywhere 
throughout the simulation domain, including the BL. His setup assumed initial net vertical field 
and is most similar to our NVF model. 

We also observe amplification of the field component in the BL in our simulations, although 
the amplification is highly time- variable. Referring back to Fig. [3j we see a bump in within the 
BL for the ZNF case (panel c), which clearly indicates shear amplification of the B^ component. 
There are also bumps in B z and B^ in panels a and b of Fig. [3j respectively. However, since panels 
a and b correspond to simulations having net vertical and azimuthal flux, respectively, the bump in 
B z and B^ within the BL is due, at least in part, to a combination of flux freezing and compression 
of accreted material from the disk on the surface of the star (j|7.1|). 



To disentangle amplification of magnetic field due to flux freezing and compression from shear 
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time 

Fig. 4. — Panels a,b,c show a:B,eff for simulations M9d,e,f respectively. Note the different scale on 
the time axis in the figures. 

amplification, we plot the maximum (in zu) of the r.m.s averaged field (in z and <j)) as a function 
of time in Fig. [5] for simulations M9d,e,f. We see that the NVF and ZNF simulations (black and 
red curves) exhibit peaks, and the NVF simulation (black curve), in particular, has one large peak 
at t ~ 560. The NAF simulation (blue curve), on the other hand, exhibits a steady increase and 
contains small local peaks. The steady increase can be understood as compression of advected 
azimuthal field within the BL. 

The peaks in Fig. [5] are due to transient shear amplification of the field in the BL, and the 
baseline level in the NVF simulations (black curve) is simply set by the value of in the saturated 
state of MRI turbulence in the disk. Fig. [6] shows the r.m.s value of as a function of radius 
from simulation M9d at times t — 400 when there is little or no amplification and t — 560, which 
corresponds to the peak in the black curve in Fig. [5| From Fig. [6] it is clear that amplification 
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2000 



time 

Fig. 5. — Plot of the maximum (in w) of the r.m.s average (in z and <j)) of as a function of 
time for simulations M9d,e,f (NVF, black curve; NAF, blue curve; ZNF, red curve). The spikes 
correspond to transient shear amplification of the component of the field within the BL. 



of the field does indeed take place in the BL and that it is time-dependent. The cause for the 
time-dependency has not been explored, but it is possible that it is related to the re-excitation of 
acoustic modes within the BL (Q. 




Fig. 6. — Plot of the r.m.s average (in z and (/)) of B^ for simulation M9d (NVF) at two different 
points in time. The dashed line corresponds to t = 400 when there is little shear amplification 
of the magnetic field in the BL and the solid line to t = 560 when shear amplification reaches a 
maximum in time. 



Although there is field amplification in the BL in our simulations, the field always remains 
strongly subthermal, consistent with the findings of Armitage (2002). Fig. [7] shows radius-time 
plots of (equation 16]) for simulations M9d,e,f, respectively. If magnetic pressure and thermal 
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pressure are comparable, then /3 ~ 1 and if thermal pressure dominates magnetic pressure, (3 ^> 1. 
In the figures, < .06 for all time and throughout the entire simulation domain, which implies 
a significantly subthermal field. 




time 

Fig. 7. — Panels a,b,c show radius time plots of from simulations M9d,e,f respectively. There is 
evidence for amplification of the magnetic field in the BL, but the field remains strongly subthermal. 



6. Acoustic Modes in the MHD context 



Before beginning a discussion of angular momentum transport by waves in the BL and the star, 
we first show that waves are indeed present in our simulations. The waves being referred to are the 
MHD analogs of the acoustic modes discussed by Belyaev et al. ( 2012a|b ) in the hydrodynamical 
case. We start by discussing how the hydrodynamical dispersion relation is modified in the presence 
of a magnetic field and then discuss simulation results. 
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6.1. Modification to the Hydrodynamical Dispersion Relation 



Belyaev et al. (2012b) found that for the isothermal BL, there were three branches of acoustic 
modes, and they found dispersion relations that matched simulation results for each of these three 
branches. However, their results were for the pure hydro case, and we comment on how they should 
be modified to account for the presence of a magnetic field in the disk. 

It is natural to expect that the major modification to the dispersion relation of the acoustic 
modes in the limit of f3 ^> 1 is the replacement of the sound speed with the magnetosonic speed. 



s 2 + v\ (23) 



= + (24) 

where va is the Alfven speed, and 7 = 1 for an isothermal equation of state. 

In reality, the situation is slightly more complicated, and the dispersion relation depends on 
the jump in magnetic field across the BL. We show this explicitly in Appendix [Aj where we derive 
an analytic dispersion relation for the 2D modes of a compressible vortex sheet ( |Miles|1958 ; Gerwin 



1968; |Belyaev fc Rafikov|2012| ) in the presence of a magnetic field that is perpendicular to the plane 



of the modes. This is not a realistic setup, but it highlights how the hydrodynamical dispersion is 
modified when a magnetic field is introduced. 

The main result of the analysis, in the context of astrophysical BLs, is that in the limit of 
^> 1, the dispersion relation for the model system is only modified by terms of O This 
result, which may be expected to hold for more general field configurations, shows that when gas 
pressure dominates magnetic pressure, acoustic modes should still be present, and their properties 
should only be slightly modified as compared to the hydro case. In a typical astrophysical context, 
such as in MRI turbulence, the magnetic field is subthermal, meaning f3 > 1. As we have shown in 
£j5j the same is also true inside the BL, meaning that the modes present in the layer and its vicinity 
are essentially acoustic modes, weakly modified by a magnetic field. 

In addition to magnetosonic waves, an ideal fluid in the presence of a magnetic field also 
supports Alfven and slow waves. One may wonder whether there are also modes of the system that 
correspond to a coupling of these waves in the disk across the BL to either gravity or sound waves 
in the star. This subject is beyond the scope of the present work. However, we mention that in our 
simulations, the modes that show up most clearly are the analogs of the acoustic modes discussed 
by Belyaev et al. ( 2012a|b ). Thus, it is natural to expect that these magnetosonic modes are the 



most important ones for transporting angular momentum in the presence of a magnetic field, when 
/3>1. 
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6.2. Simulation Results 

We observe excitation of magnetosonic modes in our simulations, which persist, typically, for 
the duration of the simulation. Fig. [8] shows both the simple vertical average and the cross-section 
through the z — plane of both w-^fpv^ and the magnetic field magnitude for simulation M9f at 
time t = 1000. The dominant mode in Fig. [8] is the m = 14 lower branch mod^J It appears 
prominently in zu-^/pv m (panels a and b), but only weakly in the magnetic field magnitude (panels 
c and d). The measured pattern speed of the m = 14 mode from simulations is ftp « .4, which 
is in good agreement with the theoretically predicted pattern speed for the m = 14 lower branch 
mode, ft P = Alb ( |Belyaev et aL]|2012b| ). 



We also point out that the mode depicted in Fig. [8] is effectively two-dimensional, since there 
is not much difference between the vertically averaged image of w-^fpv^ (panel a) and the cross- 
sectional slice (panel b). The magnetic field, on the other hand, exhibits much more vertical 
structure (panels c and d), which is indicative of turbulence. 

The lower branch has an evanescent region in the disk around the corotation radius and waves 
incident on the evanescent region are reflected back towards the BL. The edges of this evanescent 
region are located approximately at the Lindblad radii, which are implicitly defined by the condition 
^(^lr) = =t ^(^lr)M, where k is the epicyclic frequency. The dashed lines in Fig. [8] show 
the edges of the evanescent region, as predicted by hydrodynamical theory assuming a Keplerian 
rotation profile, and the solid line shows the corotation radius of the mode in the disk. Even 
though we now have MRI turbulence in the disk, shocks still reflect off inner edge of the evanescent 
region. This fact, together with the agreement between the measured pattern speed and that from 
hydrodynamical theory, means hydrodynamical theory provides a good description of the properties 
of magnetosonic modes even in the presence of MRI turbulence in the disk. As argued earlier, this 
can be attributed to the fact that < 1 in the simulations, so the magnetic field does not 
significantly affect the properties of acoustic modes. 

Although the hydrodynamical and MHD cases exhibit many similarities, some differences are 
apparent. The first difference is that MHD simulations are typically "noisier" than hydro runs 
( Belyaev et al.|2012a|b ) with a superposition of modes present. The second difference is that unlike 



the hydrodynamical case, there can be significant wave action tunneling through the evanescent 
region in MHD runs. 

The presence of a superposition of acoustic modes in the disk is a generic feature, independent 
of initial magnetic field geometry. In addition to the lower branch, which is persistent in all of our 
simulations for t > 60, we also observe the upper branch at the beginning of some simulations, but 
only in transience. Fig. [9] shows the upper branch at t = 40 in simulation M9c. The mode depicted 
in the figure has m = 14 and a measured pattern speed of flp = .84. This agrees well with the 



1 For a discussion of wave branches, including the lower branch, see 



Belyaev et al. 



(2012b) 
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Fig. 8. — Snapshots at t = 1000 from simulation M9f. Panel a) shows the simple vertical average 
of zu-y/pv^j, and b) shows a slice through the z = plane of w^/pv^. Panel c) shows the simple 
vertical average of the magnetic field magnitude and d) shows a slice through the z = plane of 
the magnetic field magnitude. The dashed lines mark the edges of the evanescent region for the 
m = 14 lower branch mode, and the solid line marks the location of the corotation radius in the 
disk for that mode. 



theoretically predicted pattern speed for this mode, which is ftp = .87. 
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Fig. 9. — Image of the simple vertical average of w^fpv^ at t = 40 for simulation M9c. The m = 14 
upper branch mode is clearly present in the image. The magnetic field is not shown, since at this 
time, the MRI instability has not yet fully developed. The dashed lines mark the edges of the 
evanescent region for the m — 14 upper mode, and the solid line shows the corotation radius for 
that mode. 



7. Angular Momentum Transport by Waves in the BL, Star, and Inner Disk 

Having shown that magnetosonic modes are present in our simulations and have similar prop- 
erties to the acoustic modes discussed in Belyaev et al. (2012b), we now show that they are able 
to transport angular momentum in the star and inner disk just as in hydrodynamical simulations. 



7.1. Density Gap in the Inner Disk and Accretion onto the Star 

We begin our discussion of wave transport of angular momentum by showing the evolution of 



the density profile. Fig. 10 shows radius time plots of Eo(tx7, t) in simulations M9d,e,f. In each 
case, a gap in the density develops in the inner disk in the vicinity of the BL {w > 1). The gap 
is opened up relatively rapidly, over a time of < 10 orbits, and in panels a and c, gap formation is 
followed at a later time by gap deepening. Concentrating on panel c, after the initial gap opening 
event at t w 200, the density profile of the gap stays approximately constant, until t « 1100 when 
it undergoes substantial gap deepening. For t > 1100, the gap is gradually filled in by the action 
of turbulent viscosity in the disk. 



Belyaev et al. (2012b) also observed gap formation in their hydrodynamical simulations, which 
was driven by angular momentum transport by acoustic modes. Consistent with their results, we 
find that magnetosonic modes are responsible for opening up the gap in our MHD simulations, and 



we discuss the link between acoustic modes and the gap in the inner disk in more detail in £7.2 



Material gets removed from the inner part of the disk between the BL and the wave reflection 
point at the inner Lindblad resonance by the action of magnetosonic modes. In this region gas orbits 
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Fig. 10. — Panels a,b,c show radius time plots of from simulations M9d,e,f (NVF, NAF, ZNF) 
respectively. A density gap in the disk near w — 1 develops in all of the simulations and in panels 
a and c gap deepening is observed at a later time. 



the star faster than the pattern of acoustic modes, and every time a fluid element encounters a weak 
shock (into which modes evolve) it loses angular momentum and moves towards the star, giving rise 
to continuing accretion. Fig. [TT] shows plots of the instantaneous mass accretion rate through radius 
vj, M(w), for simulations M9d,e,f (panels a,b,c) and for a M = 9 hydro simulation from Belyaev et 



al. ( |2012b ) (panel d). The solid curves correspond to periods when magnetosonic/acoustic modes 
are excited to a high amplitude (high shock accretion rate) and the dashed curves to periods when 
their amplitude is not as high (low shock accretion rate). Each of the solid curves in panels a-c 
corresponds to one of the gap formation or opening events in Fig. [lOj 



In panels a-c, the dashed and solid curves converge to approximately the same value of M in 
the outer disk in Fig. [TTJ but the solid curves have a significantly higher amplitude than the dashed 
ones in the inner disk [w > 1). The extra contribution to M for the solid curves in the inner disk 
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Fig. 11. — Panels a, b, c show plots of the instantaneous mass accretion rate, M, from simulations 
M9d,e,f (NVF, NAF, ZNF) at the times indicated in the legends. Panel d shows the mass accretion 
rate for an M = 9 hydro simulation considered in Belyaev et al. ( 2012b[ ) (their 3D9a simulation). 
The solid and dashed lines show M during periods when acoustic modes have relatively high and 
low amplitudes, respectively. The large spike in accretion rate in the inner disk for the solid curves 
as compared to the dashed ones is caused by shock accretion in the inner disk due excitation of 
modes to a high amplitude. 



in panels a-c is provided by magnet osonic modes, and the base level of M in the outer disk is due 
to MRI turbulence. 



Evidence for this assertion is provided by considering panel d of Fig. [TTJ Since the simulation 
corresponding to panel d is hydrodynamical, there is no MRI in the disk. This explains why the 
dashed line is everywhere approximately zero, since in the absence of waves excited to a high 
amplitude or MRI, there is no efficient mechanism for angular momentum transport in the disk. 
On the other hand, the shape and amplitude of the solid curve in the inner disk in panel d are 
similar to those of the solid curves in panels a-c. This shows that in the MHD case, the mass 
accretion rate in the inner disk due to shocks can dominate that due to MRI, during periods when 
modes are excited to a high amplitude. 



Indeed, we can write the following expression for the combined mass accretion rate close to 
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the star (interior to the evanescent region in the disk) driven by both the MRI and shocks: 



M 



d In (aqs 2 Zzu 2 ) 



dlnw 



+ 



m 



3(2 - q) 



+ (it) 



(25) 



where q(w) = —dlnQ/dlnw is the local shear rate (q = 3/2 for a Keplerian disk but is a function 
of w for a general rotation profile), m is the azimuthal wavenumber of the mode, and AE/E is the 
density contrast across the shock (mode) fronts. The first term in brackets is the contribution from 



MRI (a here is solely due to MRI) and can derived using the results of Lynden-Bell &; Pringle ( 1974) 



and |Rafikov (2012). The second term is the contribution from magnetosonic waves following from 



Belyaev et al. (2012a), namely from their equation (B7). For a characteristic azimuthal wavenumber 
of m ~ 15 — 30 and MRI viscosity a < 10 -2 (see Fig. [I]) measured in our simulations, one finds 
that both transport mechanisms provide comparable contribution to M if AE/E ~ 0.2. 

However, during periods of strong wave excitation AE/E increases and the wave term strongly 
dominates M because of the steep (cubic) dependence of M upon AE/E. For instance, the am- 
plitude of the m = 14 magnetosonic mode changes by a factor of ~ 5 in simulation M9f between 
t = 900 and t = 1100. This results in a predicted mass accretion rate due to shocks in the inner 
disk that is higher by approximately two orders of magnitude at t = 1100 as compared to t = 900. 

This effect of enhanced mass accretion rate by shocks during periods when magnetosonic modes 



are excited to a high amplitude is demonstrated by the solid curves in Fig. 11, Equation (25) also 



explains the characteristic spatial dependence of M on ro: approximate conservation of the wave 
angular momentum flux (neglecting for the moment dissipation at the shock fronts) results in (see 
equation (25) of Belyaev et al. ( 2012a| )) 



AE 



cx 



\n(w)-n P \ 



1/2 



(26) 



This expression shows that AE/E is highest near the BL, where \tt(zu) — Qp\ is large and goes to 
zero in the evanescent region where ft = ftp. This explains why during the high wave amplitude 
episodes, when mass accretion is clearly dominated by dissipation of acoustic modes, M is largest 
near the BL, but rapidly decays to small value as w approaches the evanescent region in the disk. 



Note that inside the BL MRI does not operate, a = 0, and the first term in equation (25) 



vanishes, leaving wave dissipation as the only means of mass transport. On the other hand, in 
the disk outside the corotation radius, waves do not contribute to transport and second term in 
equation (25) disappears. Thus, MRI remains the only mass transport mechanism outside the 



resonant cavity in which the pattern of standing acoustic modes exists, in agreement with the 



picture outlined in Belyaev et al. (2012a) and Belyaev et al. (2012b). 
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7.2. Angular Momentum Transport due to Waves and MRI 

In this section, we connect the density gap in the inner disk to magnetosonic modes and discuss 
angular momentum transport in the star, BL, and inner disk more generally. 



Fig. 12 shows radius time plots of the angular momentum current, C5, from simulations 
M9d,e,f. In the disk, far away from the BL, angular momentum transport is due primarily to MRI 
turbulence (Q, and Cs is approximately constant at a given radius once the MRI has saturated. 
We note that the striations appearing in Cs are caused by waves. The striations in the panels 



of Fig. 12 are predominantly downward sloping in the inner disk, indicating inward propagating 
waves. These waves are excited at the "front" of MRI turbulence which propagates radially outward 
through the disk over the course of the simulation. Although potentially a nuisance, these waves 
do not play a major role in angular momentum transport in our simulations. 

In contrast to the outer disk, where the value of Cs at a given radius is approximately constant 
in time once MRI turbulence has saturated, Cs varies by more than an order of magnitude in 
the star and inner disk over the course of a simulation. This implies a mechanism of angular 
momentum transport that is decoupled from MRI turbulence. As we shall shortly demonstrate, 
this mechanism is transport of angular momentum by magnetosonic modes. Given this fact, we 
can directly link magnetosonic modes to the gap formation and deepening events in Fig. [lOj since 



these are simultaneous with the periods when Cs is large in the BL and star in Fig. 12 



In addition to waves dominating angular momentum transport in the star and BL, we find 



that they also contribute significantly to angular momentum transport in the inner disk. Fig. 13 
shows Cs and superimposed on the same plot for simulations M9d,e,f at times t = 550, 250, and 
1100, when the angular momentum current due to waves is high in each simulation. The dashed 
vertical line denotes the radius at which Cs = 0, and for an anomalous viscosity model, Cs = 



and dVtjdw — occur at the same radius (equation [16]). However, we see in Fig. 13 that the 
radius at which Cs — is significantly displaced into the disk compared to the radius at which 
dVtjdw — 0. This implies that Cs contains a contribution due to waves in the disk, since Cs < 
for outgoing waves up to the corotation radius in the disk, which is located at zu cr w 1.8 (m = 14 



mode) in each of the panels in Fig. 13 



We conclude that when waves are excited to a high amplitude, they significantly modify the 
profile of Cs in the inner disk compared to what is expected for a turbulent viscosity model. In 
particular, if both magnetosonic modes and anomalous viscosity contribute appreciably to Cs, then 
the radius at which Cs — in the disk is neither at the corotation radius of the mode, nor at the 
radius where dVtjdw — 0, but somewhere in between. 

We also mention that the variation of Cs with radius in the star and BL in each of the panels 
in Fig. 13 is consistent with Cs(w) for the lower branch (Fig. 12 of Belyaev et al. ( 2012b[ )). Since 



the lower branch is directly observed in our simulations (q6j), this provides further confirmation 
that angular momentum is transported by waves in the BL and star. 
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Fig. 12. — Panels a,b,c show radius time plots of the angular momentum current Cs from simula- 
tions M9d,e,f (NVF, NAF, ZNF) respectively. Notice that Cs changes sign in the vicinity of the 
BL with Cs < for w < 1 and Cs > for w > 1. 



7.3. Hydro dynamical Nature of Angular Momentum Transport in the Star and BL 



Having established that magnetosonic modes transport angular momentum in the BL and star, 
and also influence transport in the inner disk, we now demonstrate that the mechanism of angular 
transport by acoustic modes is essentially hydrodynamical in nature. This provides evidence for the 
argument presented in £6T that when (3 ^> 1, the properties of acoustic modes are not significantly 
affected by the magnetic field, and the results of Belyaev et al. (2012b) are valid. 
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(c) 

Fig. 13. — Panels a,b,c show plots of Cs and for simulations M9d,e,f at times t — 550, 250, and 
1100. The dashed vertical line denotes the radius at which Cs — and is significantly displaced 
into the disk relative to the radius at which d£l/dzu = 0. In general, when angular momentum in 
the disk is transported by both a magnetosonic mode and turbulent viscosity, the location at which 
Cs = occurs between the radius at which dVtjdw — and the corotation radius of the mode, 
vj ct « 1.8 in this example. 



To show that stresses in the BL are hydrodynamical, we split Cs into magnetic and nonmag- 
netic components (equation fl4 



Fig. 



14 



shows Cs,b an d Cs,h for each of simulations M9d,e,f, 
and several things are immediately apparent. The first is that C$,b is several times larger than 
Cs,h m the disk, which is consistent with the fact that the magnetic stress is several times the 
hydrodynamical stress for MRI turbulence ( Sorathia et al.||2012 ). The second is that Cs,b vanishes 
going into the star and there is no bump in Cs,b within the BL. This means that magnetic stresses 



are not important for transport in the BL. This is consistent with the results of Pessah & Chan 



(2012), who found that in the shearing sheet approximation, the magnetic stress oscillates about 
zero for a rising rotation profile {dVt/dw > 0), as in the BL. The third is that \Cs : h\ ^> \Cs,b\ in 
the star and BL, which is definitive proof that a hydrodynamical mechanism of angular momentum 
transport operates there - namely acoustic modes weakly modified by a magnetic field. 
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Fig. 14. — Panels a,c,e show radius time plots of Cs,b from simulations M9d,e,f (NVF, NAF, ZNF) 
respectively, and panels b,d,f show radius time plots of Cs,h from simulations M9d,e,f. The scale, 



colormap, and axes labels for each simulation (M9d,e,f) are the same as in Fig. 12 Notice that 
Cs,b vanishes inside the star, and there is no bump in Cs,b within the BL. On the other hand, 
Cs,h, is large inside the BL and is both large and approximately constant (spatially) within the 
star, which suggests transport by waves. 



7.4. Boundary Layer Width 



Belyaev et al. (2012a) and Belyaev et al. (2012b) found that in hydro simulations the BL 

s 2 /g(vcj*)). Using their definition for 



settled to a radial width of 6 — 7 radial scale heights (h. 
the BL as the region in which 

0.1 < (v^(vj))/v k {vj) < 0.9, 



(27) 



we confirm that after the initial gap opening event, the BL in our MHD simulations also settles 
down to a radial width of 6-7 radial scale heights (Fig. [Ts] ) . However, the BL undergoes widening 
events that are coincident with gap deepening events (j |7.1| ) and the periods when Cs is large in the 
star and BL (j |7.2| ). Stochastic widening events were also observed by |Belyaev fc Rafikov (2012) in 
2D hydro simulations, and it appears that in both hydro and MHD cases these events are mediated 
by the action of shocks on the fluid in the inner disk. 



The fact that the BL widening events occur simultaneously with gap deepening events, suggests 
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a common origin, implying that magnetosonic modes play a critical role in regulating the width of 
the BL. However, given the occurrence of stochastic widening events in the simulations, the width 
of the BL is time-variable and it is not possible to assign a single number to it. Nevertheless, what 
the simulations suggest is that the BL initially settles down to a steady state width of 6-7 radial 
scale heights and is subsequently widened, during periods when magnetosonic modes are excited 
to a high amplitude. 
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Fig. 15. — Boundary layer width in units of scale height s~ 2 (defined implicitly in equation 
as a function of time, for each of our simulations. The black, blue, and red curves correspond to 
NVF, NAF, and ZNF simulations respectively. The solid lines correspond to simulations with the 
higher resolution for a given magnetic field geometry. 



8. Discussion 



We have shown that even in the presence of MRI turbulence in the disk, waves dominate 
angular momentum transport in the BL, and influence transport in the inner disk. The importance 



of waves for angular momentum transport in the BL was previously shown by Belyaev et al. (2012b) 
for the purely hydrodynamical case. 



We find that when a magnetic field is introduced, the acoustic modes discussed by Be lyaev 



et al. (2012b) become magnetosonic modes. These magnetosonic modes have a dispersion relation 
that differs from the dispersion relation for acoustic modes by factors that are Since 
the field in our simulations is highly subthermal < .06) everywhere, including the BL, the 



hydrodynamical theory of Belyaev et al. (2012b) remains a good approximation to the dynamics 



of magnetosonic modes in our simulations. 

Unlike MRI turbulence, which can be parametrized by an anomalous viscosity, the transport 
of energy and angular momentum in the BL is a non-local process, as has been emphasized in 



Belyaev et al. (2012b). Energy and momentum are carried away from the BL by acoustic waves 
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weakly modified by magnetic field into both the star and the disk and can be deposited far away 
from the BL. This has important implications for stellar spinup, heating of the star, disk and the 



BL itself, and for observational manifestations of the BL phenomenon (Belyaev et al. 2012b). In 



the MHD case, additional dissipation is produced in the disk due to the operation of MRI but, as 
one can infer from Fig. [I] where olb is seen to vanish in the BL, the MRI-generated heat is unlikely 
to penetrate into the BL and affect its structure. 

Amplification of the magnetic field in the BL that we see in our simulations (§51) is consistent 



with earlier theoretical expectations (Pringle 1989) but does not result in efficient magnetically- 



mediated transport of angular momentum in this region. Our study is in agreement with the 



analytical work of iPessah & Chan (2012), who analyzed the evolution of MHD modes in presence 



of shear with dVt/dw > and demonstrated that although the energy density of these modes can 
be amplified significantly, their associated stresses oscillate around zero, rendering them inefficient 
at transporting angular momentum. As a consequence, angular momentum transport in the BL is 
mediated by quasi-acoustic modes, and in particular shock dissipation of these modes in the inner 
disk. 



Previously, Steinacker & Papaloizou (2002), Armitage (2002) performed isothermal MHD sim- 



ulations of the BL, and the setup of Armitage (2002) in particular was nearly identical to ours. 



Therefore, a natural question that arises is why did these authors not observe acoustic modes and 
recognize their importance for angular momentum transport in the BL? Interestingly, the answer 
to this question may be that the modes were in fact observed in these early simulations, but their 
presence and significance were not fully appreciated. 



For instance, Fig. 3 of Armitage (2002) shows an image of relative density perturbation aver- 



aged over the z-dimension. The density fluctuations in the image exhibit a leading/trailing mode 
pattern that looks strikingly similar to the lower mode trapped between the BL and the evanescent 
region in the disk, see Fig. [8j Note that AS/S oc E -1 / 2 for an acoustic mode propagating into the 
star by angular momentum conservation, see equation (26), which explains why the amplitude of 



density fluctuations rapidly decreases going into the star in Fig. 3 of Armitage (2002). 



It is harder to comment on whether Steinacker & Papaloizou (2002) observed acoustic modes 



since they used a different initial setup than ours and prescribed a rotation profile for the BL region 
in their simulations. In contrast, the profile of in the BL is naturally established in our simulations 
due to action of the sonic instability early on in the simulations (Belyaev et al. 2012a|b ), which 
widens the initial interfacial region into a self-consistent BL. Nevertheless, it is quite possible that 



Steinacker & Papaloizou (2002) also observed acoustic modes excited in the BL, as they briefly 



mentioned the presence of "stochastic spiral patterns" in their simulations, which could well have 



been the manifestation of acoustic modes (e.g. see Fig. 7 of Steinacker & Papaloizou (2002)). 



Another important question that arises is given that we have used an isothermal equation of 



state, what might be expected for a more realistic equation of state? As discussed in |Belyaev et 



al. (2012b) one pathology of the isothermal equation of state is that the Brunt- Vaisala frequency, 
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N 2 = (7 — l)g 2 /s 2 , is equal to zero. This means that our computational model does not support 
gravity waves and in a real star with nonzero Brunt- Vaisala frequency, magnetosonic waves in 
the disk could couple to gravity waves in the star increasing the number of possible modes in the 
star-BL-disk system. 

Another difference is that the sound speed is not constant in a real star and is an inverse 
function of the radius. This means that the cutoff frequency below which sound waves do not 
propagate is also an inverse function of the radius and without a damping process, such as radiative 
damping, waves propagating into the star would reflect back towards the BL. A simple way to 
model this effect in simulations would be to introduce a reflecting wall at the inner boundary of the 



simulation domain. Glatzel (1988), Belyaev & Rafikov (2012) found analytically that for the case 



of a supersonic shear layer with a linear velocity profile, instability was present for both reflecting 
wall and radiation boundary conditions. Therefore, although reflection at an inner boundary would 
likely modify the details of angular momentum transport by waves, in particular the spinup rate 
of the star and possibly the wavelength of the dominant mode, magnetosonic modes would still be 
excited in such a system. 
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A. Dispersion Relation for the Compressible Vortex Sheet in the Presence of a 
Constant Magnetic Field that is Oriented Perpendicular to the Motion 

Here we present the dispersion relation for modes of the compressible vortex sheet with 
wavevectors lying in the x-y plane, when a z-directed magnetic field is present. Under this set 
of assumptions, the B • VB term in the MHD equations is zero, which means Alfven waves are 
not permitted, the slow wave is null, and only the fast magnetosonic wave survives. In fact, in 
such a setup the magnetic field simply provides pressure support and behaves as though it were 
a fluid with an effective adiabatic index of 7^ = 2. Although this is not a realistic description 
for astrophysical boundary layers, it gives us a flavor for how the dispersion relation of acoustic 
modes is modified when a weak magnetic field is introduced. In particular, it shows that for a weak 
magnetic field, the corrections to the dispersion relation are small and are of order O 

Since all fluid quantities in our setup are independent of z we are essentially studying a 2D 
problem. In that case, the ideal MHD equations with an adiabatic equation of state can be formu- 
lated as 

^ + V-(pv) = 0, (Al) 
d{pv) + V.(pvv) = -Vfp+^V (A2) 



dt vr J V 

P = Kp\ (A3) 

B z = ap, (A4) 



where, we have replaced the induction equation ^ with the frozen-in law ( A4) and introduced the 



"magnetization" of the fluid per unit mass, a. Equations ( Al||A4 ) can be simplified to 

g + V-(pt;) = 0, (A5) 
d ( pV Kv-(pvv) = -v(k P ^ + ^p 2 ). (A6) 



dt vr 7 V 2/i 

Written in this form, it is clear that for our setup the magnetic field behaves in the same manner 
as a fluid with 7^ = 2. 



We now describe the compressible vortex sheet setup, following the notation of |Belyaev fc 



Raflkov (2012). The velocity profile of the model setup is given by 



W-{-C*<°o, (A7) 

where V y is a constant, and a velocity discontinuity is present at x = 0. For x > 0, the density, mag- 
netization, and sound speed are constant and are given by p + , cr+, and s+, respectively. Similarly, 
for x < these parameters are constant as well and are given by p_, cr_, and s_, respectively. 
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Using the techniques laid down in ( Gerwin|1968 ; Alexakis et al.|2002 ; Belyaev fc Rafikov|2012 ), 
it is a straightforward if lengthy exercise to derive the dispersion relation for linearized perturbations 
of this system. Thus, we omit the derivation and simply state the dispersion relation, which is 



e 2 [(M ms + ^) 2 e" 1 (l + S)- 1 - 1] (M ms - ^>) A = [(M ms - if) 2 - 1] (M ms + ip) 4 



(A8) 



Here, e = p+/p~ is the density contrast across the vortex sheet, tp = u/k y s+ is the dimensionless 
phase speed, M ms = V y / s \ + v / 
the relation 



is the magnetosonic Mach number, and 5 is defined through 



e(l + S) 



As in the text, va denotes the Alfven velocity. 



s - + v X- 



(A9) 



Equation (A8) is identical to equation (39) of Belyaev & Rafikov (2012) for the unmagnetized 



case, except that the Mach number has been replaced by the magnetosonic Mach number, and 
there is an additional factor of (1 + The first of these modifications, the replacement of the 

Mach number with the magnetosonic Mach number, is to be expected and was argued for in the 
text. The second is harder to anticipate and a nonzero 5 can be thought of as coming from a 
difference in the effective 7 of the fluid across the vortex sheet. The effective 7 of the fluid is a 
suitably weighted combination of 7 and 7^. 



The major point of equation (A8) as regards astrophysical BLs is to show that for a weak 



magnetic field, the dispersion relation for the compressible vortex sheet is only modified up to 
terms of O (v\/s 2 ) (i.e. O Since the dispersion relation for the three wave branches of 

acoustic modes discussed in ^6] is simply related to the dispersion relation for the compressible 
vortex sheet ( |Belyaev et al.|[2012a ), we reason that properties of acoustic modes are also only 
modified by terms of O in the presence of a magnetic field. Since (3 tends to be large in our 

simulations, even in the BL, we reason that the magnetic field due to MRI perturbs the acoustic 
modes, but does not modify them in a fundamental way. 



